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THEORY AND RESULTS ON COLLECTIVE AND COLLISIONAL EFFECTS 


FOR A ONE -DIMENSIONAL SELF-GRAVITATING SYSTEM 

By Frank Hohl 
Langley Research Center 

SUMMARY 

The equilibrium properties of one- dimensional self- gravitating systems are inves- 
tigated analytically. One-dimensional models are used to perform computer experiments 
tracing the evolution of stellar systems. The stationary solution of the Vlasov equation 
for a one- dimensional system of stars as computed for an interesting class of initial con- 
ditions is found to correspond to a minimum -energy configuration. The results of the 
numerical experiments are compared with theory. For initial energies far from the 
minimum equilibrium energy the system becomes unstable and breaks up into smaller 
clusters. A variational principle was applied to the one- dimensional stellar system to 
show that stationary distribution functions which decrease monotonically in going outward 
from the center of the system are stable. That other stationary distributions may be 
unstable is illustrated by means of computer experiments. The one- dimensional model 
is of interest as an approximation to the distribution of velocity and mass normal to the 
galactic plane of a greatly flattened galactic system. Observational results for the gravi- 
tational force normal to the galactic plane of the Galaxy agree with the results obtained 
from the one-dimensional model. Thermalization effects for systems containing small 
numbers of stars were investigated. The fluctuation of the kinetic energy was found to be 
inversely proportional to the square root of the number of particles in the system. 

INTRODUCTION 

Many problems in stellar dynamics involve phenomena occurring in inhomogeneous 
systems in which the interaction between the particles is fully described by a self- 
consistent field operating in phase space. Because the particles interact by means of the 
long-range Coulomb force, each particle is under the simultaneous influence of a large 
number of other particles. Therefore, stellar systems will respond to any perturbation 
in a collective manner, and a study of such systems is concerned essentially with the 
N-body problem. 

The collective phenomena do not depend on two-body collisions such as occur in 
ordinary gases, and therefore the collective effects will be present in collisionless 


systems. Since the number of particles in the system is large, a distribution function 
can be used to describe the density of particles in phase space. The distribution function 
must then satisfy the Vlasov equation (the self-consistent set of the Maxwell equations 
plus the collisionless Boltzmann equation). In using the Vlasov equation to describe a 
stellar system the number of masses which make up the system is assumed to become 
infinite while the total mass remains constant. Although such an approach allows 
description of the system by means of a distribution function which must satisfy the 
Vlasov equation, solutions to the time-dependent nonlinear Vlasov equation are, in gen- 
eral, very difficult to obtain. An attempt is therefore made to condense the 10 H stars 
which a galaxy may contain into about 10^ superparticles. Numerical or computer models 
are then used to perform computer experiments simulating Vlasov phenomena by fol- 
lowing the simultaneous motion of a number of superparticles of the order of 10^. One- 
dimensional models are used to effectively solve the time- dependent nonlinear Vlasov 
equation for various systems. 

Computer models used to study plasmas and stellar systems are essentially of two 
types. The first type is the Lagrangian model in which the self-consistent motion of 
large numbers of particles in phase space is followed. In such a model the number of 
particles which can be treated is limited by computer storage and time. The second type 
of computer model is the Eulerian model in which the system behaves like a fluid in 
phase space. The macroscopic quantities describing the system are then obtained by 
solving the appropriate partial differential equations on the computer. In the present 
study the main interest is in systems which can be described by the Vlasov equation. 
However, it was found advantageous to use a Lagrangian model to simulate the time 
development of the systems investigated. Thus, the main interest is not in solving the 
Vlasov equation for the macroscopic quantities describing a system but rather in per- 
forming numerical experiments by solving the N-body problem. 

Direct solutions of the N-body problem have become possible only with the advent 
of high-speed electronic computers. Since the pioneering work of Pasta and Ulam 
(ref. 1), many workers have used computer models in their investigations of the N-body 
problem. 

The simplest system of stars which can be considered is a one-dimensional system 
with the stars stratified into plane parallel layers. All parameters then vary only in 
one direction. With the exception of the work done by Lecar (ref. 2) with small numbers 
of sheets, the author is not aware of any previous numerical experiments to simulate 
stellar systems with a one- dimensional sheet model. Henon (ref. 3) used a computer 
model of concentric spherical mass shells to study the dynamical mixing of spherical 
star clusters. Recently, Hockney (ref. 4) and Hohl (ref. 5) have used two-dimensional 
rod models to simulate a cylindrical galaxy. The results obtained with the two- 
dimensional model show that filamentary and spiral structure can be obtained by purely 
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gravitational effects. However, these calculations should be repeated by using point 
masses confined to move in a plane. Lindblad (ref. 6) calculated the two-dimensional 
motion of a number (up to 192) of mutually attracting mass points in a given central fielcL 
of force. By placing the mass points initially in a system of concentric rings with circu- 
lar velocities, Lindblad investigated the mutual disturbances in such a system to simulate 
the spiral structure of galaxies. The remaining numerical calculations in stellar dynam- 
ics appear to be three-dimensional calculations for systems with small numbers of stars. 
Thus, von Hoerner (refs. 7 and 8) made extensive calculations with up to 25 stars. In 
such a system the effects of the individual encounters are large and cannot be smoothed 
by averaging over many calculations. Aarseth (ref. 9) performed similar three- 
dimensional calculations with up to 100 stars and applied the results to clusters of gal- 
axies. The irreversibility in stellar systems with up to 32 stars was investigated by 
Miller (ref. 10). 
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SYMBOLS 

constant defined by equations (19) and (42) 
acceleration 

matrix elements defined by equation (85) 

energy correlation function 

speed of light in vacuum 

Debye length, v^,(47rGp) 

gravitational field 

energy distribution function 

distribution function 

gravitational constant 

total energy density 

potential energy density 
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S'p 

H 

M 


m 

N 


n 

P 


r 

T 

t 

U 

v + ,v_ 

V n 

v 


kinetic energy density 
Heaviside unit step function 
mass 

mass per unit area 

total number of particles in system 

particle density 

potential energy 

radial coordinate 

kinetic energy 

time 

1 9 

total particle energy, -mv^ + mcp 

contours defined in figure 2 
normalized equilibrium contour 

velocity 

thermal velocity 


W = T + P 


x 

X’ 

z 

a 


position coordinate 
position defined by cp(x') = e 0 
dummy variable 

constant defined as 167rGxm^A\ L 

V2m k 


and variable used in equation (66) 


4 
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$ 

<p 

X 

X' 

<> 


Dirac delta function 

potential energy defined by equation (22) 
energies defined by equation (118) 
integration variable 


function used in variational method 
«x 
-x c 


variable, ^ V + (x)dx 


thermal energy 

relative fluctuation of kinetic energy 
energy ratio 
mass density 
dimensionless time 
characteristic period, (4 ttGp) 


1/2 


dimensionless gravitational potential 


gravitational potential 


dimensionless position coordinate 


dimensionless position coordinate defined by <3?(x') 


time- averaged quantity 


Subscripts: 


eq 


equilibrium 




summation indices 


min 


minimum 


o initial 

s system boundary 

Superscript: 

k summation index defining waterbag contours or ordering index 

THEORY OF SELF-GRAVITATING STELLAR SYSTEMS 

The analytical portion of the present study is confined to an investigation of steady 
states for stellar systems. A computer model is used to simulate the time behavior of 
the system and to study its approach to an equilibrium state. The simplest star system 
is one in which the stars are stratified into plane parallel layers. The velocity distribu- 
tion and density then vary only in one direction. For such a model the stars can be repre- 
sented by a large number of mass sheets. The motion of such sheets is sufficiently 

simple that the trajectories of several thousand stars can be followed on an electronic 

computer. Thus numerical experiments to investigate the evolution of a stellar system 
can be performed and the results compared with the equilibrium conditions obtained from 
the Vlasov equation. 

As indicated by Michie (ref. 11) the importance of orbital (or phase) mixing in the 
initial evolution of a system of stars needs to be investigated. As the system of stars 
evolves, the gravitational field changes with time and the stars follow complicated tra- 
jectories along which the individual stellar energies are not conserved. Henon (ref. 3) 
has recently performed numerical experiments with a system of concentric spherical 
shells to study the relaxation of the mean gravitational field and the resulting approach 
to equilibrium for a spherically symmetric star cluster. The analogy between models in 
which the material is stratified in parallel planes and those in which the material is dis- 
tributed in concentric spherical shells has been examined by Woolley (ref. 12). 

According to Oort (ref. 13) the velocities of stars normal to the galactic plane are 
decoupled from the other velocity components. Since the force on a star is approximately 
normal to the galactic plane, the stars will oscillate perpendicular to the galactic plane 
with a period independent of the revolution of the galactic system. A one-dimensional 
model representing a system stratified in infinite parallel planes can therefore be used 
as an approximation to the distribution of velocity and density of stars normal to the 
galactic plane of a greatly flattened galactic system. Camm (ref. 14) has considered 
steady solutions of the collisionless Boltzmann and the Poisson equations for such a model, 
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and his results compared favorably with the observed densities. More recently Lindblad_ 
(ref. 15), Woolley (ref. 16), and Oort (ref. 17) have used a one- dimensional analysis to 
study greatly flattened galactic systems. Prendergast (ref. 18) has investigated the gen- 
eral solutions of the steady- state one- dimensional self-gravitating system. 

The present study is concerned primarily with the conditions under which a system 
will approach a stationary state and for which a system will become unstable and break 
up into smaller systems. The relaxation, or characteristic, time for a system of stars 
is given by 

r c = (AirGp)' 1 ' 2 (1) 

where G is the gravitational constant and p is the mass density. The characteristic 
length of interest is the Jeans or Debye length D which is defined by 

D = v T r c (2) 

where v T is the thermal velocity of the system of stars. The dimension of a system of 
stars near equilibrium is of the order of a Debye length. Thus, for a three-dimensional 

O 

system the number of stars within the Debye sphere is very large, that is, nD 0 » 1 
where n is the average density of stars. For example, a galaxy contains somewhere 
from 10® to 1C)12 stars so that nD® ~ 10® to 1C)12. Similarly, for the one -dimensional 
model the number of stars in a Debye length is large, or nD » 1. The effects of colli- 
sions between individual stars can then be neglected, inasmuch as the distance traveled by 
a star before collisions significantly modified the star's trajectory is much larger than 
the dimension of the system (that is, the Debye length of the system). Thus, an actual 
stellar system as well as the one-dimensional model is well within the Vlasov regime. 

For nD » 1, the numerical experiments using the one- dimensional model should simu- 
late the time development of the distribution function as given by the Vlasov equation. 

It should be mentioned that computer calculations such as those performed by von 
Hoerner (refs. 7 and 8) and by Aarseth (ref. 9) involving the three-dimensional motion of 
up to 100 stars are not described by the Vlasov equation. 

In the one-dimensional model the forces between the sheets are long range and the 
sheet model includes individual as well as collective behavior of the stars since the com- 
puter simply solves the equations of motion of the stars in the system. The exact behav- 
ior of a system depends on the graininess and will be affected by going to the "fluid limit" 
as implied by the Vlasov equation. As shown in a subsequent section, equivalent models 
with different graininess (or collisions) give identical results for any time scales of pres- 
ent interest; thus, graininess has a negligible effect on the system. Also, since the col- 
lective behavior is not affected in going to the fluid limit, good agreement can be expected 
between the "experimental" results and the theoretical results obtained from the Vlasov 
equation for the steady state of the system. 
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The phrase "metaequilibrium" refers to the equilibrium state which is established 
by the interaction of the stars with the smoothed potential of the system in a time of the 
order of 2 ttt c . The metaequilibrium can be considered a steady state only as long as 
graininess or collisions can be neglected. On a long time scale, binary encounters 
eventually cause the system to approach a state of statistical or thermal equilibrium. 

An actual three-dimensional stellar system can never completely attain a state of statis- 
tical equilibrium since this would involve a Gaussian velocity distribution. The absence 
of a potential barrier outside an actual stellar system permits the escape from the sys- 
tem of all stars with positive energies. Thus, the escape of stars will prevent the estab- 
lishment of a Gaussian distribution necessary for a state of statistical equilibrium. 
Nevertheless, the time constant with which the three-dimensional system will tend toward 
a state of statistical equilibrium (even though it can never be completely attained) is of 
the order nD^r c , where n is the density of stars. This time is of the same order of 
magnitude as the time constant involved in Chandrasekhar's (ch. n of ref. 19) calculation 
of the dynamical friction. For the one-dimensional model the behavior is different 
because of the potential barrier outside the system. Fluctuations due to graininess may 

now cause the establishment of a state of statistical equilibrium in a time which is at 

o 

least of the order of (nD) t c . 

Lecar (ref. 2) has performed numerical integrations of a one- dimensional system 
of stars and has shown that no thermalization exists to order nDr c . For a plasma with- 
out an external electric field this has been shown analytically by Eldridge and Feix 
(ref. 20). Presently, only times of order r c are of interest, during which collisional 
effects will be negligible. 

Comparison of the gravitational fields acting on stars in a one-, two-, and three- 
dimensional system shows that they are very similar with the exception of the field acting 
on stars near the boundary of the system. Figure 1 shows the gravitational field for a 
one-, two-, and three-dimensional system with constant density. Because of the l/r^ 
dependence of the field outside the three-dimensional system, stars can escape; this is 
not possible for the one- dimensional system. However, inside the system the force on a 
star in all three systems is proportional to the distance from the star to the center of the 
system. 


Stationary Solutions 

In problems concerning the structure and evolution of stellar systems Camm 
(ref. 21) has shown that in the Vlasov limit (nD » 1), the difference in mass and struc- 
ture of individual stars can be neglected. Therefore, all the star masses are set equal 
to m so that a distribution function f(x,v,t) completely defines the state of the system. 
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When nD » 1 the effects of graininess are negligible and the 
fies the usual equation of stellar dynamics, that is, the Vlasov equation 


9f 9f „ 8f . 

~ + v — + E-h = 0 

9t 9x 9v 


0(j9 

where E = - is given by the Poisson equation 
ox 


8 % = 4irGm 
9x 2 



( 3 ) 


(4) 


The limitations of the Vlasov equation in describing stellar systems are discussed by 
Kurth (ref. 22). 

In a study of the time development of a system of stars the distribution function at 
t = 0, that is, f(x,v,t=0), could be used as the initial condition in an attempt to solve the 
time-dependent nonlinear Vlasov equation. But this problem has not yet been solved 
analytically for general initial conditions. Therefore, a high-speed computer is used to 
determine the evolution of the system from Newton's laws instead, and equations (3) 
and (4) are used only to obtain the steady- state solution for the system. 


The Virial Theorem 


The virial theorem can be applied to the one -dimensional system to obtain a rela- 
tion between the potential and kinetic energy of the system in equilibrium. Consider an 
arbitrary volume of phase space which is convected with the flow. The integral over this 


volume ^ x^f(x,v)dx dv is a well-defined function of time, and according to the trans- 
port theorem (p. 131 of ref. 23) 


-~2 §§ x 2f (x,v)dx 


dv = ^ ^g[x^f(x,v)]dx dv 


(5) 


13 3 8 9 j-vf 

where the substantial derivative — = -rr + v — + a — . By the Liouville theorem ^ = 0, 

Dt 9t 9x 9v Dt 

and thus equation (5) when multiplied by m/2 becomes 


t $ II A <*> v > dx dv = f II '^^ 2 )* 1 dv 


= ^ mv^f(x,v)dx dv + mxaf(x,v)dx dv 


( 6 ) 
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If the region of integration spans the entire system of particles, the integral 
^ x%(x,v) dx dv is constant and the left-hand side of equation (6) vanishes. Also 

^ mv%(x,v)dx dv = 2T (7) 

where T is the total kinetic energy of the system. Using p(x) = \ mf dv and a = - 4^ 

J dx 

m equation (6) gives the expression 

0 = 2T 


- 


dx 


The Poisson equation can be written as 


P(x) = 4^tf 

Therefore, an integration by parts yields the result 


( 8 ) 


(9) 




4irG -Xg dx 2 dx 


X 

8ttG\ 


-X c 


. j- r Xs 

8ttG J_ x ^dxj 


( 10 ) 


The last expression is simply the potential energy P of the system, normalized so that 
P = 0 when all masses are at the same point. Equation (5) then becomes 

2T - P = 0 (11) 

irrespective of the distribution function. In most of these numerical studies it was found 
that the system very quickly approached conditions such that equation (11) was satisfied. 


In the steady state, 



Waterbag Distribution 

and the Vlasov equation takes the form 


v — + E — = v — - ^2 — 
9x 9v 8x dx dv 


= 0 


( 12 ) 


Using the method of characteristics to solve the partial differential equation given by 
equation (12) results in the subsidiary equation 

dx _ dv 
v dcp/dx 


(13) 
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which gives the result 

i mv^ + mcp = U = Constant (14) 

where U is the total energy of a star. Therefore, any solution of equation (12) haS7.thaez^ri-_i 
form 

f (x, v) = F(U) (15) 

Thus, if a time- independent equilibrium state exists 

f(x,v,t-°°) = F(U) (16) 

where F(U) may be any function of U. Of course, only functions F(U) that are stable 
are of interest in this study. In general, the form of F(U) depends on the initial dis- 
tribution and must be obtained by following the time development of the Vlasov equation. 
However, there is one type of initial distribution for which F(U) is known without 
actually solving the Vlasov equation. For this distribution the initial f is taken to be 
constant over a certain region of phase space and is zero outside this region. Figure 2 
illustrates the distribution. According to equation (12) f remains constant along the 
different trajectories so that the region can only change its shape with time while keeping 
its area constant. For this reason the distribution function just described has been called 
the waterbag model by DePackh (ref. 24). When DePackh considered the waterbag model 
he was interested primarily in a solution for the linearized oscillations about the equilib- 
rium state. However, a more interesting application of the waterbag model seems to be 
in the present study of the nonlinear problem. The waterbag model is of interest in this 
study primarily because it allows calculation of the exact equilibrium configuration of 
the one-dimensional star gas for comparison with the computer results. 

Two initial distributions considered in the present study are the waterbag model 
and a distribution which has a constant density over a region of the X-axis and has a 
Maxwellian velocity distribution. For the waterbag model, the stationary distribution 
function F(U) is known to be constant for 0 ^ U = me and to be zero for U outside of 
this region; me is some maximum energy to be determined later. For the second dis- 
tribution, F(U) is assumed to be Maxwellian. 

The equilibrium solution for the waterbag model is obtained as follows. The initial 
shape of the waterbag for most of the calculations is taken to be a rectangle defined by 
the area in phase space between ±x Q and ±v Q . The function F(U) is now used in the 
Poisson equation (eq. (4)) and the integration over dv is changed to an integration over 
dU. Since the density at a given x only is of interest, dU = mv dv or 

JTT 

dv = uu is obtained. The Poisson equation then becomes 
\|2m(U - mcp) 
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( 17 ) 


d 2 y _ BTOmA f™ dU _ A ^ 

dx 3 V2m VLT - m cp 


where me is the energy such that F(U) =0 for U > me and F(U) = A for 
0 ^ U ^ me . Since the area of the system in phase space remains constant, the value of 
A for the initially rectangular waterbag is obtained from the relation 


or 



dx dv = N 


(18) 


A = 


N 

4x 0 v 0 


(19) 


A first integration of equation (17) gives the result 

v2 


'dc£^ _ _ 8 \/2 7 rGNm 
dx 


3x o v o 


c?) 3 /2 _ e 3/2 


f 


( 20 ) 


where <p = ^ = 0 at x = 0 has been chosen to determine the constant of integration, 
dx 

The condition ( =0 is simply the statement that the force on a star vanishes in a 

\dx/ x= o 

plane dividing the system into two equal masses. A second integration gives the final 
result 


±x = 


1/2 

3 Vo \ r%3/2 _ (e _ e) 3/2j- 


.8 \J2 7rGNmy 


XT- 


1 / 2 . 


d? 


( 21 ) 


Let x s be the coordinate defining the boundary of the one- dimensional system; 
then <p(x s ) = e and = 27rGNm. If these values for <p and d^/dx are used 

in equation (20), the result is 

1 2/3 


e = 


f37rGNmx 0 v 0 . 


2^2 


( 22 ) 


The value of x s is obtained from equation (21) by taking e as the upper limit of the 
integral. Thus 


/ 2 2\ 1/3 

The value of the maximum velocity is given by the expression 

v s = \/27 = ^37rGNmx 0 v 0 ^ 1//3 


f n 2 2 

9x*v* 
o o 

V 7rGNm , 


, 1/3 


(23) 


( 24 ) 
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( 25 ) 


Equation (17) shows that the particle density is given by 

n(x) = -=*— (e - <p) l/2 

V2 x 0 v 0 

where use was made of equation (19). 

If the dimensionless variables 

X = (27rGNm/e)x (26) 

and 

$ = £ (27) 

are used in equation (21), the simplified equation obtained is 



The resulting <f> and E as functions of y are shown in figure 3. 


Multiple- Contour Waterbag Distributions 

Another waterbag distribution which was investigated numerically is one with a hole 
in the center, referred to as a two- contour waterbag. The only possible F(U) is then 
one which has a constant value A for me 0 g U g me where me is the maximum and 
me 0 is the minimum star energy in the system and F(U) is zero outside this range. 

The Poisson equation then takes the form 

87rGm A $ (^e - <p - ]je 0 - cp^j (x < x') (29) 


and 


= 87rGmA\l2 
cbt 

where x' is given by <p(x') = e 0 . A first integration gives the result 

^2 


(x > x') 


and 


(sf) 2 = " x^ GmA [ (e " <p)3/2 + ° 2 ] (x > x,) 


(30) 


(h^) = " (e " <p)3/2 “ ( e ° ' <p ) 3/2 + Cl (x < x,) (31) 


(32) 
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= 27rGNm 


Using the boundary conditions cp(0) = 0, 0 = °’ ^( Xs ) = e ’ and 

yields the constants Cj and C2 

c , - -( e 3/2 - e „ 3/2 ) 


c 2 = 

the final result 


' £ 3/ 2 . So 3/2) 


second integration gives 

MVf *GmA ) 1/2 J 0 ^ 3/2 - (e ~ C) 3/2 ] -[e 0 3/2 - (e 0 - ?) 3/2 ]} 


(33) 


(34) 


±x= 32 

l 3 


1/2 


d? (x < 


and 


Using the 


and 


- 1/2 

?|I T G mA )' 1/2 J[ e 3/2 . < £ . c) 3/2] . eo 3/2] dc (x : 


dimensionless variables 


X 


= (l |I . ca ) 1 ' 


$ 


= £ 


simplifies equations (35) and (36) to 

= lottf ' (1 ' c)3/2 ] ■ t" 372 ■ {v - c)3/2 ]} dt <X < X') 

-1/2 


±x 


and 


; 3 / 2 j d? 


** = x ’ + " (1 ‘ 

shows the normalized (dimensionless) equilibrium 
jntours are obtained from the equations 

/■4\ l~~l ‘ 


(x > x' 


e Q 

where u = — . Figure 4 shows the normalize 
four values of v. The contours are obtained 

£ } (x) = ]j2[v - <&(xjj 


and 


= ]]2[v - $(x)J 
V< 2 )(x) = ^[l - *(x)] 
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Maxwellian Distribution 


For a system at thermal equilibrium the distribution function F(U) has the form — 

F(U) = A exp(-fcU) (42) 

Following the same method used in obtaining the solution for the waterbag model yiel ds—" 
the equilibrium density 

n = 2 f f dv = 2 f A exp(-KU) - dU = 2A\I~— exp(-Kmcp) (43) 

J 0 J m cp \|2m(U - m <p) » 2mK 

Equation (4) then becomes 

2 I 

- 8irGmAlL-£— exp(-Kmcp) = 0 (44) 

^<4 iZniK 


or 


£-$-f-exp(-*) = 0 
dx z 2 


(45) 


where $ = Kxncp and a = 167rG/cm 2 A ^ 2 ^ — . A first integration of equation (44) gives the 


result 


\dxj 
m a 

x = 0. A second integration gives 


(sr) " “C 1 - exp( - 4) ] = 0 


(46) 


d4> 


where the constant of integration was determined by the boundary condition — = 0 at 


±x = — In 
\Ja 


1 + ^1 - exp(-<fr) 


1 - \J 1 - exp(-4>) 


Solving equation (47) for gives the equation 


4>(x) = -In 


1 - tanh 


2/fa 



One of the boundary conditions to be satisfied is 

.. d<i 
lim — 

dx 
or 


lim = 27 rGNm 2 K 


- 2 «r 


a = (27rGNm‘ 

The solution is then given by the equations 

<z>(x) = - In sech(7rGNm 2 /<x) 
^ ' xm ' ' 


(47) 


(48) 

(49) 

(50) 

(51) 
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and 


E(x) = -27rGNm tanh(7rGNm 2 /<x) 


(52) 


n(x) = ^ 7 rGN 2 m 2 K sech 2 (7rGNm 2 Kx) (53) 

Ll 

where 1/ k is the average kinetic energy per star at any given position. This particular 
solution has also been considered by Camm (ref. 14). 

Figure 5 shows the variation of 4> and E for the Maxwellian distribution where 
the values 47 tG =1, m = 1, k = 1, and N = 4 were arbitrarily chosen. In figure 6, 
the variation of the observed gravitational field normal to the galactic plane of the Galaxy 
is compared with the results obtained from equations (20) and (52). These values of 
E(x) normal to the galactic plane were obtained by the astronomer Oort and have been 
tabulated by Lindblad (ref. 15). The variation of E(x) for the Maxwellian distribution 
was found to be in good agreement with the experimentally obtained variation. 


Energy Considerations 

The steady- state solution for the waterbag distribution was obtained by using con- 
servation of area in the two-dimensional phase space. For the waterbag model it was 
possible to obtain F(U) corresponding to a given initial condition. It must therefore be 
determined whether energy considerations allow the relaxation from a function f(x,v,t=0) 
of two independent variables to a function F(U) of the energy only. Thus, if an equilib- 
rium state is to be reached, the initial and final energy of the system must be equal. For 
the initial rectangular waterbag the kinetic energy T is given by 

T(t=0) = Jf i mv 2 f(x,v,t=0)dx dv = f ° C ° v 2 dv dx = ^ v 0 2 (54) 

JJ i ° x O V 0 -Xq "V 0 ® 

As before, ±v Q and ±x 0 define the rectangular area of the initial distrioution. By use 
of equation (10), the potential energy of the initial waterbag is found to be 


x 2 

P(t=0) = — (27TGNm) 2 - — C °( 27rGNm £-\ dx = f 7rGN 2 m 2 x n (55) 
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The total initial energy W is then given by 


W(t=0) = T(t=0) + P(t=0) = fl(|v 0 2 + 27rGNmx 0 J 


(56) 


For a given constant value of A = 
is easily found to be given by 


N 


4x 0 v 0 


, the minimum value that equation (56) can attain 
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( 57 ) 


The total energy of the equilibrium state that is defined by equation (2 1) is now 
computed. The kinetic energy is given by 


-eq 


2 N f x S C r 

mv f dx dv = - \ \ yU - mcp dU dx 
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where dv = 


dU 


. = was used and, as before, <i> = cp/e and y = (27rGNm/e)x. By 

v2m(U - mcp) 

substituting equation (20) into equation (10), the potential energy for the equilibrium water- 
bag becomes 


p e „ = 3L(2*GNm) 2 - JL f Xs (i|lGNm\r e 3/2 _ 
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The equilibrium state then has a total energy given by the expression 


^eq ^eq + ^eq “ g 
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dx 


where 


f Xs (l - $) 3 / 2 dy -0.988 


(60) 


For arbitrary x Q and v Q it is clear that equations (56) and (60) cannot be equal. 
The equilibrium energy W e q depends only on the product x Q v 0 whereas W(t=0) is 
the sum of two independent quantities that are functions of the two independent variables 
v Q and x Q . The total energy as well as the area in phase space must be conserved. 
Thus, the only state compatible with the conservation of area in phase space and with the 
dynamics of the Vlasov equation does not have the same energy as the initial state. 
Therefore, a steady state cannot exist for an initially rectangular waterbag distribution. 
However, the equilibrium state has the interesting property that it is a minimum- energy 
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configuration. As is shown in the section on computer simulation, the system will do its 
best within the limitations of energy conservation to reach the equilibrium state. 

The Minimum- Energy Principle 

The minimum- energy property of the equilibrium waterbag distribution can be 
demonstrated analytically. First consider the initial rectangular waterbag distribution. 
Equation (57) gives the minimum initial energy for such a distribution. The ratio of this 
minimum initial energy to the equilibrium energy is 


W 


mm 


(t=0) 


w, 


eq 


1.03 


and the rectangular waterbag with the minimum energy has still more energy than does 
the equilibrium waterbag. The minimum- energy property can be generalized to any 
shape of the waterbag distribution by showing that the equilibrium state has the least pos- 
sible energy. 

Consider the waterbag distribution shown in figure 2. Let points along the contour 
be described by V + (x) and V_(x) where the plus sign indicates the upper contour and 
the minus sign indicates the lower contour. The equations of motion for V ± (x) are 

8V. 8V. 

+ r- E (61) 


For the equilibrium state 


~atT 


0 and equation (61) becomes 


dV. 


*f- E = 0 


(62) 


A variational method is used to demonstrate the minimum- energy property of a one- 
dimensional self- gravitating system. This method has previously been applied by Hohl 
and Staton (ref. 25) to a plasma described by a single -contour waterbag distribution. 

The use of the variational method requires the energy density g in the system. Using 
the first integral in equation (10), the potential energy per unit length gp is given by 


gp = -2mAxV + 



(63) 


where 6 = f V + (x)dx 

-*s 

A 3 

g T is simply j mV + . 
has the form 


and 4irGm^ - 2A0j = E. The kinetic energy per unit length 
Thus, the energy density for the present problem g(x,V + , 
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g = j mV + 3 - 2mAxV + pirGm^ - 2A6>) 


L 

where A is the magnitude of f . Thus, the total energy of the system is given by 


w = j* S g(x,V + ,e)dx 

-X c 


(64) 


(65) 


For the sake of simplicity the contour of the waterbag has been taken to be sym- 
metric so that V + = ~V_. Following Courant and Hilbert (ref. 26), the extremum of 
equation (63) is found by the usual methods of variational calculus. The end points ±x s 
of the contour are held fixed and the function V + (x) receives a variation coj(x) where 
a is an arbitrary constant and r?(x) is an arbitrary function which vanishes at the end 
points. It is then found that 


9W 

3a 




d a g 3g' 
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where 


r x s 

\ 77 dx = 0 

J -x s 


(67) 


Since 77 is an arbitrary function satisfying equation (67), the condition for obtaining an 
extremum is 


d / 9g\ 3g 
dx\3 V + j 30 


= 0 


( 68 ) 


If the expression for g given by equation (64) is used in equation (67), the result is 

dV, 

v+ "df‘ E '° (69) 

which is identical to equation (62) for the upper contour. Application of the Legendre 
criterion of the second variation of g can tell something about the extremum just cal- 
culated. For the problem under consideration 

a 

^-4 = 2mAV + (70) 

3V“ 

and, since for small perturbation about the equilibrium state V + is positive (or zero), 
the Legendre criterion shows that the extremum cannot be a maximum. Thus, the mini- 
mization of the total energy of a waterbag distribution leads to the contours satisfying the 
stationary Vlasov equation. The consequence of the minimum- energy property is that 
starting from any nonequilibrium state, energy conservation will prevent the system from 
reaching the steady state described by equations (21) and (25). This result was to be 
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expected from the analysis of DePackh (ref. 24) which shows that, for a plasma, small per- 
turbations of the equilibrium state are not damped. Nevertheless, the interesting point is 
that numerical results show that the system does its best within the limitation of energy 
conservation to approach the equilibrium state. In general, it is found that the equilibrium 
state is approached very closely whenever the initial energy is not too different from the 
energy of the equilibrium state given by equation (60). 


Generalized Minimum-Energy Property 

The minimum- energy principle can be extended to arbitrary distribution functions. 
The waterbag model illustrated in figure 7 is used in the analysis. The contours V^(x,t) 
and V_ y (x,t) bound surfaces of constant f = f k . According to equation (3) the phase 
fluid bounded by the contours is incompressible. In the limit of a very large number of 
contours the waterbag model can be used to approximate arbitrary distribution functions. 

The distribution function describing a multiple- contour waterbag distribution has the 

form 


f = ^ A k h(v - V^) - h(v - V^) 


(71) 


where H(z) is the Heaviside unit step function. This distribution function satisfies the 
Vlasov equation so that the contours are streamlines of the flow in phase space, that is 
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(72) 


The equilibrium contours are given by 

(k) 8V i k) 

vi K; -TP- - E = 0 (73) 

± dx 

The force per unit mass E is obtained from the equation 


HI = -4irGm J f dv = -4irGm ^A k fv| k ^ - V^j 
cu J k ' ' 


(74) 


where it is assumed that there is a total of N stars, each of mass m, in the system. 

fkl AA (kl 

To simplify the equations, symmetric contours V_^ ' = -VI ' = v' ' are assumed in the 
derivation of the minimum- energy property. The results obtained are not affected by 
such an assumption. 


A new variable is now defined as 
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(75) 
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where x);' is the end point of the contour k. In terms of 0 y ' the gravitational field 
s 

is given by 


E(x) = 47rGm 




(76) 


Let g be a function such that the total energy W of the system is given by 
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g(x, 6>( k ) jV ( k ))dx 


(77) 


Taking the extremum of the integral for W, subject to the condition that 

N =I 2 v (k) (4 k) 

k v 

is a constant, requires that the contours v^ satisfy the Euler- Lagrange equations 


(78) 


0g* d 8g* 


90 (k) dx 9 y(k) 


= 0 


(79) 
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where g* = g + 2XA k v v ' and X is, as yet, an tmdetermined Lagrangian multiplier. 

Since the end points are held fixed, the condition on the total number of particles is 
redundant and X = 0. If the end points x^ are not held fixed in the variational process, 

o 

additional equations will appear which will not change the main results given. 

For a multiple waterbag the total kinetic energy per unit length is 




A k v ^ 


k k 

From equation (10) the potential energy per unit length is 
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Therefore, the expression for g is 
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The Euler-Lagrange equations then become 
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dx 
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= 0 


r (k) dv 
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dx 


E = 0 


(83) 


(84) 


Equations (84) are the equations for the equilibrium contours satisfying the Vlasov 
equation. 

If equations (84) are to represent a minimum- energy configuration, then the Legendre 
criterion of the second variation of g must be satisfied. Following Courant and 
Hilbert (ref. 26), the Legendre criterion for the case of several unknown functions 0^ 
is that the quadratic form whose coefficient matrix has the elements 


— 


9 2 g 


13 0v(^3v(i) 

must not be negative. For the present problem, only the diagonal elements of equa- 
tion (85) are nonzero and are given by 

,00 


(85) 


a kk = 2mA k v ' 


( 86 ) 


Since 2mv 


(k) . 


is never negative, the Legendre criterion requires that 

A k >0 


(87) 


From the definition of f given by equation (71) it can be seen that equation (87) is equiv- 
alent to stating that the distribution function must always decrease in going outward from 
the center of the system where f = f-^ must be the largest. If equation (87) is satisfied, 
the system is a minimum- energy configuration and is always stable. However, if Ajj > 0 
is not satisfied for all k, the system is not a minimum- energy configuration. If the 
Legendre criterion is not satisfied, the system may be unstable since two or more con- 
tours can now be deformed while keeping the total energy of the system constant. In a 
subsequent section the results of a numerical experiment are presented which show how 
two- contour systems with Aj = -A2 become unstable. 


COMPUTER SIMULATION OF STELLAR SYSTEMS 


The simplest system of stars that can be considered is the one- dimensional system 
in which the stars are stratified into plane parallel layers and all parameters vary only 
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in one direction. The nonrelativistic motion of stars in one dimension is simple enoughs 
so that the motion of large numbers of stars can be followed on an electronic computer. 

Initial distributions of the waterbag type were used extensively in the numerical 
experiments since they allow a comparison with the results obtained in the preceding 
sections. Two types of initial waterbag distributions have to be considered separately. 
The first type is a "fat" waterbag for which the dimension of the system is near D and 
the total energy is near the energy of the equilibrium distribution. The system then 
quickly approaches the equilibrium state with the exception of the arms in phase space 
which develop to accommodate the excess energy. The second type of initial distribution 
is a "thin" waterbag which has a total energy much larger than the energy of the equilib- 
rium distribution. For this second type the system cannot even come near its equilibrium 
state and an instability develops which causes the system to break up into small clusters. 


Description of the One-Dimensional Model 


The one -dimensional model consists of a system of N stars which are represented 
by mass sheets. These sheets are of infinite extent along the Y- and Z-axes and the 
sheets are constrained to move along the X-axis. The equations of motion of all the N 
stars are solved simultaneously by computing the gravitational field at the position of 
each star and then integrating the equations of motion for each star over a small time 
interval At and repeating the process. When two sheets meet they are allowed to pass 
freely through each other. 

The evolution of the system of stars is studied by using various nonequilibrium ini- 
tial distributions as input to the computer and then following the time development of the 
system on the computer. 

The equation of motion for the ith sheet with a mass m per unit area is given by 
the expression 

A (t) 

-52- - E (v) < 88 > 


where Xj is the position of the ith sheet and E(x,t) is obtained from 

N 

- * E 8 M = 47rGp(x,t) = 47rGm ^ 6jx - Xj(t)j 

j=l 

The expression for the gravitational field is then 


E(x,t) = 27rGm 
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(89) 
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where 


sgn(z) = -1 
= 0 
= 1 


(z > 0) 
(z * 0) 
(z < 0) 


Figure 8 illustrates the variation of the gravitational field for an eight-sheet system. In 
the actual numerical calculations the field was computed by first ordering the mass sheets 
according to their coordinate so that 


x (h) g x (k+l) 


(91) 


The gravitational field is then given by 


(92) 


and the summation indicated in equation (90) does not have to be performed at every time 
step and for each star. The sorting routine used in ordering the stars is very fast and 
takes advantage of the fact that advancing the motion of the stars for a small time inter- 
val will not get the stars too far out of order. 


The motion of the mass sheets can also be computed by a more exact method. That 
is, the times at which neighboring pairs of sheets cross are computed and the shortest of 
these crossing times is used to recompute the new positions and velocities for the sheets 
affected. The time for the next crossing is then found and the process is repeated. 

The accelerations are constant until two sheets cross. This very accurate program 
has been used by Lecar (ref. 2) to study certain "invariants" of the system. The exact 
program is also used in the present study to investigate thermalization for systems with 
small numbers (up to 40) of stars. However, the exact program takes a large amount of 
machine time and is not suitable for investigating the motion of large numbers of sheets. 


The electronic data processing system was used to calculate the self-consistent 
motion of systems containing several thousand stars. The position and velocity of each 
sheet are computed at successive times tj, t 2 , t 3 , . . ., tj. For each time step 
At = t n+ i - t n the new positions and velocities of the sheets are computed from the 
equations 
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where 


d2 Xj( tn ) _ 
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-A-J = 27rGm 2 H[xi(t n ) - Xj (t n )J 


d x i( tn ) 1 d x i( t n) d x i( t n- 1) 

dt 3 At|_ dt 2 dt 2 J { } 

3 / 3 

Introduction of the term d Xj/dt allows use of a larger time interval At in the com- 
putations while keeping the error very small. The accuracy of the computer program 
was checked by comparing the results for various values of At and also by checking 
the reversibility of the motion. For example, the error in the computed total energy of 
the star system was less than 0.005 percent for the systems investigated. Similarly, the 
total momentum of the system remained constant to within 0.0005 percent. If n is the 
average density of stars and D is the Debye length or dimension of the system, then the 
numerical experiments using the model described are very nearly exact for nD » 1 and 
actually simulate the solution of the nonlinear Vlasov equations. 

The potential energy for a system of N stars, each of mass m per unit area, is 
given by the equation 

r n 

P ‘ 8iG < 2lrGNm > 2 (*N - x l) - 2 - *j-l) < 97 > 

j= 2 

where Ej(x) is defined by equation (90) so that for each j, Xj_ j < x < Xj. This defini- 
tion of the potential energy is such that if all stars are at the same point, then the poten- 
tial energy of the system is zero. The kinetic energy of the system is given by the usual 
expression 

N 

T = | m J vj 5 (98) 

j=l 

The total energy Uj of a star is of interest and is given by 

u j = | mv? + m(?j (99) 

The gravitational potential <p ^ is obtained after first ordering the mass sheets 
according to their increasing x- coordinates as indicated by equation (91). Thus, 
indicates the gravitational potential of the sheet labeled j (which is the first sheet) 
Dropping the subscript j results in 
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where N is an even integer. Also, cp was chosen equal to zero at the midpoint of the 
system which divides it into two equal masses. 


The virial theorem can also be applied to a discrete system in equilibrium. Equa- 
tion (6) then takes the form 
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where 
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= 0 and m^) = m for all k. Use of the expression 
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yields 


^ mx^ d ^ = 27rGm 2 [(N - l)x^ + (N - 3)x® + (N - 5)x^ + . . . 
k=l dt 
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If symmetry of the system is assumed such that x' 
be written in the form 
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, equation (105) can 
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The potential energy of the system is given by 
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Comparison of equation (107) and equation (106) shows that 

Y mx^ k ) d2>: ^ - -P (108) 

k=l d ‘ 2 

and therefore 

2T = P (109) 

as was shown before for a continuous system. 

Initial Distributions Near Equilibrium 

Since the computer simply solves the equations of motion of the stars in the system, 
the one-dimensional sheet model includes individual or star- star interactions as well as 
the interaction of the stars with the smoothed potential of the system. The graininess or 
collisional effects are affected by going to the "fluid limit" as implied by the Vlasov 
equation. 

The effects of graininess on the time development of a system can be determined by 
varying the number of sheets while keeping the total mass of the system constant. The 
effects of graininess must be checked to determine whether the model is adequately 
described by the Vlasov equation. Figure 9 shows the time development of two equivalent 
systems with the same initial distribution but with different graininess. The curves indi- 
cate the variation in time of the kinetic energy for an initially rectangular waterbag 
distribution. The initial conditions are obtained by using a random number generator to 
give a nearly uniform distribution over a rectangular region in phase space. For all the 
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numerical calculations the gravitational constant was normalized such that 4nG = 1. 
Also, all time scales are normalized to r c 0 = r c (t=0). The initial ratio of kinetic to 
potential energy for both curves in figure 9 is 2, whereas the equilibrium value is 0.5. 
The upper curve is for a system of 1000 stars with nD ~ 1000 and m = 2. The lower 
curve is for a system of 2000 stars with nD ~ 2000 and m = 1. The oscillations of the 
kinetic energy are identical for the two systems, an indication that the graininess effects 
do not affect the development of the system. 

Comparison of the curves in figure 9 shows that after only a few periods nonlinear 
damping or phase mixing has almost eliminated the oscillations in the kinetic energy and 
the kinetic energy remains near the equilibrium value. These results indicate that the 
Vlasov equation correctly describes the system. 

To determine when the effects of graininess become important, the evolution of sev- 
eral equivalent systems with varying numbers of stars is compared. Figure 10 shows 
the variation of the kinetic energy for equivalent systems which contain from 20 to 
500 stars. For the system with 500 stars the behavior is still very similar to that shown 
in figure 9. The time behavior of the energy distribution function and the density for 
equivalent systems which contain more than 500 stars is also very similar. As can be 
seen in figure 10 , when the number of stars of the system becomes 200 or less the time 
behavior of the kinetic energy differs markedly from that of figure 9. Also, the time 
development of the energy distribution function and of the density is different for the sys- 
tems shown in figure 10. In general, it was found that oscillations in the various param- 
eters of the system persisted much longer when the number of stars of the system was 
less than 500. 

Figure 11 shows the time development of the density of the waterbag distribution for 
a system of 2000 stars with a mass m per unit area equal to one. The ratio of initial 

energy to equilibrium energy f W(t— 0)\ j or S y S ^ em j S 133 and the initial ratio of 

V W eq / 

kinetic to potential energy is 2. The corresponding variation of the kinetic energy is 
shown by the lower curve of figure 9. The dashed curve in figure 11 is the theoretical 
equilibrium density as obtained from equation (25). After only a few periods 2ttt c , 0 , 
the "experimental" density is close to the calculated value. The time development of the 
energy distribution function F(U) for the same system is shown in figure 12, where 
F(U) represents the relative number of sheets in an energy interval AU. The dashed 
line shown in figure 12 represents the theoretical distribution. Again, after only about 
six periods the experimental distribution closely approximates the theoretical distribu- 
tion. Thus, even for values of 1.33 for the ratio of initial to equilibrium energy the sys- 
tem approaches its nonstationary equilibrium state closely in a time of the order of 
2ttt c 0 . The best way to see the approach to equilibrium is to take a sequence of photo- 
graphs showing the time development of the system in phase space. Figure 13 shows 
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such a sequence for the system of 2000 stars. Each star is represented by a small__ 
circle. While approaching the nonstationary equilibrium state indicated by the area _j 
inside the oval at time t = 33.0r c 0 in figure 13, the system rotates in phase space with 
a frequency near 1/27 tt c . An interesting feature of the time development shown in fig- 
ure 13 is the development of "spiral arms" in phase space. Such arms were found to 
develop in all systems having an initial energy larger than the equilibrium energy. The 
relatively few particles in the long arms shown in figure 13 are required to accommodate 
the excess energy, and they correspond to the high-energy bump shown in figure 12 at 
time t = 38r c q . The arms develop because the period of oscillation of the sheets 
increases with energy. Figure 14 shows the period as a function of energy for the equi- 
librium state. 

The calculations just described for the 2000-sheet system with an initial waterbag 
distribution were also performed for an equivalent 1000- sheet system. The variation of 
the kinetic energy for this system is given by the upper curve in figure 9. The variation 
of the density with time is shown in figure 15 and is found to be very similar to that of 
the 2000- sheet system (fig. 11). The corresponding time development of the energy dis- 
tribution function shown in figure 16 is also very similar to that of the 2000- sheet sys- 
tem (fig. 12). 

The results of Henon (ref. 3) show that the oscillations for a system of spherically 
concentric mass shells are damped more rapidly than those for the one- dimensional sys- 
tem shown in figure 9. The reason for this is that the period of oscillation of a spherical 
shell increases much more rapidly with increasing energy than it does for a mass sheet 
in the one -dimensional system. Therefore, the effectiveness of phase-mixing, which 
depends on the change in period with changing energy, is different for the two models. 

In the equilibrium state the gravitational potential is time independent and the tra- 
jectories of individual stars form fixed closed contours in phase space as shown in fig- 
ure 17. Figure 18 shows the actual trajectories for five mass sheets for the time- 
dependent system shown in figure 13. The phase- space trajectories shown in figure 18 
are typical for systems which are initially near the equilibrium state. Initially, while 
the system is approaching the equilibrium configuration the orbits are very much per- 
turbed; this causes the redistribution of the orbits required to approach equilibrium. 

After only about two orbits the trajectories are nearly the same as the theoretical orbits. 

The time behavior described so far is characteristic of the waterbag model whenever 
the initial energy given by equation (56) is near the equilibrium energy given by equa- 
tion (60). Since the system approaches the equilibrium configuration within the limits of 
energy conservation, it appears that the steady state given by equation (21) has some sig- 
nificance. The stability of the steady state has been checked by use of an initial distribu- 
tion which satisfies equation (21). Even after many periods 2irr c , the system is found 




to remain in the equilibrium state. The effect of a perturbation of the equilibrium state 
has also been checked. Figure 19 shows the time development of a local sinusoidal per- 
turbation of the equilibrium state. Because the outermost stars are now outside the main 
system they have a longer period and cause the development of arms. After this small 
perturbation the nonstationary equilibrium state can never be completely reached. 

A system with an initially constant density between x 0 and -x Q and with a 
Maxwellian velocity distribution has a time development similar to that for the waterbag 
distribution. For example, if the initial ratio of kinetic to potential energy is of the order 
of the equilibrium value, then the system shows a time behavior very similar to that 
shown in figure 13. However, the arms in phase space develop more rapidly because 
some of the high-velocity mass sheets are now able to stay outside the main system 
longer and therefore have a longer period. For example, figure 20 shows the variation 
of the kinetic energy of such a system with an initially Maxwellian velocity distribution. 
The initial ratio of kinetic to potential energy is 2. As shown in figure 20, the kinetic 
energy of the system quickly approaches its equilibrium value. The corresponding time 
development of the density and energy distribution function is shown in figure 21. The 
dashed curves are those obtained from equation (53) and from the relation 
F(U) = A exp(-/cU). The experimental density approaches its theoretical value closely 
whereas the energy distribution again develops a high-energy bump. Figure 22 shows the 
time development of the system in phase space. With the exception of the more rapid 
development of the spiral arms the time development in phase space is close to that 
shown in figure 13 for the waterbag model. Thus, the results obtained with the waterbag 
distribution do indicate the general behavior of one-dimensional self- gravitating systems. 


Initial Distributions Far From Equilibrium 

If the energy of the initial rectangular waterbag distribution is much larger than the 
equilibrium energy given by equation (60), the system is unstable and will break up into 
smaller clusters. Figure 23 shows the time development of the kinetic energy for a 
waterbag distribution. The system contains 1000 stars and the ratio of initial to equilib- 
rium energy is 6 whereas the ratio of initial kinetic to potential energy is 32. Figure 23 
shows that the oscillations in the kinetic energy are damped only very slowly with time. 
The time development of the density for the same system is shown in figure 24. From 
the density plots it can be seen that the system quickly breaks up into three clusters 
which later combine into two clusters. Similar information is obtained from the time 
development of the velocity distribution function shown in figure 25. There are large 
fluctuations in. time of the energy distribution function of the system as shown in figure 26. 
Again, the time development of the system can best be observed by viewing it in phase 
space. Figure 27 shows the time development for the system of 1000 stars in phase 

The breakup into three clusters is clearly demonstrated. The instability which 
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develops and causes the system to break up into smaller clusters is a manifestation_of_±he- 
Jeans instability (pp. 345-347 of ref. 27). After a long time the system shown in figurei27 
condenses into two clusters with a halo of high-energy stars. In figure 28 the position ofr 
every 20th star for this system is plotted as a function of time. After only one complete 
oscillation the bunching of the stars into four distinct regions is clearly visible. 

The initial development of the instability depends on the initial fluctuation of the star 
density in phase space. To illustrate the effects of initial fluctuation of the phase- space 
density the calculations for the unstable system of 1000 stars have been repeated by using 
a different random number sequence to obtain the initial distribution. The results are 
shown in figures 29 to 32. Figures 30 to 32 show that the system now breaks up into four 
clusters. However, after a long time the system in phase space appears very similar to 
the previous case in that only two clusters remain with a halo of high-energy stars. 

Changing the graininess of the system does not affect the behavior just discussed. 

To show this, the calculations were again repeated for an equivalent system with 
3000 stars. The results are shown in figures 33 to 35. The system appears to break up 
into three clusters but it quickly condenses into two clusters by effectively winding itself 
into two spirals. The Vlasov character of the system is clearly demonstrated by the con- 
tinuity of the spiral arms shown in figure 35. 

Distributions with large initial ratios of kinetic to potential energy and which have a 
constant density between x Q and -x 0 and a Maxwellian velocity distribution are also 
unstable and will break up into smaller clusters. However, the breakup now does not 
occur as readily as for the waterbag distribution. 

Initially "Cold" Systems 

Considered next is an initial distribution in which the particles initially are placed a 
unit distance apart and have zero velocity. Each star is then attracted toward the cen- 
ter of the system with a constant acceleration. In the absence of any perturbation the 
crossings will not take place until the stars reach the center of the system. All the stars 
will reach the center simultaneously and cross at the same instant. After this time the 
force on a particular star changes sign and the system begins to expand until it reaches 
its original size with zero kinetic energy and the process is repeated. The results of the 
calculations for a 1000- star system are shown in figures 36 and 37. Figure 36 shows the 
variation of the kinetic energy with time. The curve is very similar to that shown in fig- 
ure 29 except that now there is no damping of the oscillation. If the calculation were 
extended for a longer time, the error in the numerical solution would eventually cause 
stars to get out of order and the oscillations would be damped. The variation of the den- 
sity with time is shown in figure 37. As all particles approach the center of the system 
simultaneously, the density at that point increases without limit. The system has the 



interesting property that during the expanding phase the velocity of recession of the stars 
as viewed from any particular star increases linearly with distance in either direction. 

Another initially "cold" system investigated is the two-stream case. The particles 
are initially uniformly spaced along x with half the particles at velocity v Q and the 
other half at velocity -v 0 . The variation of the kinetic energy for such a system with 
1000 stars is shown in figure 38. After a slight initial rise the oscillations of the kinetic 
energy show a strong damping. The time development in phase space for this system is 
shown in figure 39. The two streams wind around each other and in this manner attempt 
to reach some sort of equilibrium state. Figure 40 shows the evolution in phase space of 
a system which is very similar to that shown in figure 39 except that the initial positions 
of the mass sheets were random. It can be seen that the local condensations for this sys- 
tem are much more pronounced. 

In the last two systems considered, the relaxation which occurs is of a violent nature 
and it appears that such systems are suitable to check the theory of violent relaxation 
as given by Lynden-Bell (ref. 28). For the simple case of a single- contour waterbag the 
averaged distribution function <f > is given by (ref. 28) 

exp["-/3(U - /zj] 

<f > = A *=— ■ p — — 

1 + exp[-j3(U - n)j 

where U is the total energy of a star and 0 and \i are two constants determined by 
conservation of energy and mass. This distribution is formally identical to the Fermi- 
Dirac distribution. To determine whether a system approaches the distribution proposed 
by Lynden-Bell, the phase plane is divided into a number of cells and the number of parti- 
cles in each cell is found. This process is repeated for different times and the average 
of a number of such samplings gives <f>. Figure 41 shows the results for the system 
shown in figure 39. However, <f> was determined after the system had been advanced 
to t = 250 t c 0 . According to the Lynden-Bell distribution a plot of In < ^ ^ - as a 
function of U should give a linear relationship since 

-0(u - 4) - in 

These results are partially confirmed in figure 6 where the points corresponding to 

about 75 percent of the total mass lie near the dashed line. Similar results have been 
obtained by Henon (ref. 29) who used a model with concentric spherical shells to simulate 
stellar systems. 

In almost all cases considered so far, the system tried to approach a steady state 
and the dimensions of the system are then of the order of a Jeans or Debye length. This 
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result is also found if the dimensions of a globular cluster or even a galaxy are con- 
sidered. Thus, if this simple Newtonian picture is extended to a larger scale, it could 
be speculated that the present radius of, say, the universe is 


R ~ D = 


v T r c 


v t 




( 110 ) 


If a spherical geometry and uniform mass density are assumed, R can be written 


as 


R ~ 


]/gm/r 3 


( 111 ) 


or 


GM _ 1 
Rv| 


( 112 ) 


where M is the total mass in the universe. For v^ near the velocity of light c 
equation (112) becomes the well-known relation 


GM 

Rc 2 


1 


(113) 


often used (refs. 30 and 31) to define the radius of the universe. Note that as soon as 
the average density of a self-gravitating system is known, its characteristic frequency 
can be computed. For a typical globular cluster r c is of the order of 3 X 10® years; 
for a galaxy t c is of the order of 10® years. Assuming that the average density in the 
neighborhood of the Galaxy represents that of the universe, the characteristic period for 
the universe is found to be of the order of 8 x 10^® years. The characteristic period for 
the universe is close to the period of oscillation generally quoted in the oscillating- 
universe theories. 


Unstable Stationary States 

In a previous section it was shown that if a stationary distribution function always 
decreases in going outward from the center of the system, the stationary state is a 
minimum-energy state and is stable. However, if this condition is not satisfied, deforma- 
tions of the waterbag contours are possible and there is a possibility that the system is 
unstable. The interchange instability which destroys the stationary state has been 
numerically investigated for two- contour systems with Ai = -A 2 (refer to fig. 4). The 
results are shown in figures 42 and 43. Figure 42 illustrates the interchange instability 
for a system of 2000 stars with v, the ratio of the minimum to the maximum star energy 
of the stationary state, equal to 0.4. At the initial time the contours satisfy equa- 
tions (41a) and (41b). The equilibrium contours are quickly distorted while the heavier 
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outer waterbag tries to displace the inner bag. Figure 43 shows the time development 
of an unstable 2000-star system for v = 0.25. The growth of the instability is now much 
slower because the central waterbag or hole is smaller. 

Thermalization Effects In a One- Dimensional Stellar System 

In a one -dimensional system crossings will always take place so that the system can 
be expected to approach a Maxwellian distribution. An exact double-precision computer 
program was used to investigate thermalization effects for low numbers of stars. The 
exact program is one in which the shortest crossing time for a mass sheet is determined 
and the system is then advanced by that time and the process is repeated. This model, 
of course, differs from the model used for the results presented so far in that the system 
was always advanced for a constant time interval At, irrespective of the number of 
crossings during this time interval. The exact program is very accurate. For example, 
after time -reversing the numerical integration of a typical system investigated, it was 
found that the initial conditions were reproduced accurately to 12 digits. 

Lecar (ref. 2) has investigated the exact one-dimensional motion of low numbers of 

stars. He found that no thermalization occurs to order nDr c where nD ~ N. In the 

present study of thermalization effects a low number of mass sheets was chosen so that 

the system can be followed for a long period of time. In order to obtain meaningful 

results, time averages of the systems are investigated. The constant time interval used 

in the averaging process is taken to be smaller than the average crossing time for sheets 

2 2 

in the system. The time averages are taken over times of the order n D r c . It should 
be noted that the investigation of thermalization effects is still in progress and only some 
of the initial numerical results are presented here. 

For all of the systems investigated, it was found that the time- averaged potential and 
kinetic energies satisfied the virial theorem; that is, <p > = 0.5 with an accuracy of 
at least three digits. The <> signify time-averaged quantities. As discussed by Ford 
(ref. 32), the Poincare recurrence time and the general behavior of the system are depen- 
dent on the initial conditions. For example, special initial conditions can be chosen so 
that the time behavior of the system is very nearly periodic. However, for the results 
presented the initial conditions were arbitrarily chosen. 

Figure 44 shows the time- averaged velocity distribution and density for a three- 
particle system. The solid line for the velocity distribution corresponds to a Maxwellian 
distribution and is obtained by integrating A exp(-/cU) over x; that is, 
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The solid line for the density is given by equation (53). The value of 
from 

I: - 3N _ N 
2<T> + 2<P> 2<T> 


k is obtained 


( 115 ) 


The circles in figure 44 represent the time-averaged numerical results. The experi- 
mental velocity distribution is near the Maxwellian distribution. However, the variation 
of the density indicates that there is some near-periodic behavior of the system. Such 
near-periodic behavior is likely to occur for very low numbers of stars. 


For example, for a three-particle system, equation (88) can be integrated with the 
result 

x l(t) = x^O) + Vj(t=0)t + 47rGm jJ* jsgnjx^t) - x 2 (t)J + sgnJxjOO - Xg(t) j dt d? (116a) 


2 (t) = x 2 (t=0) + v 2 (t=0)t + 4?rGm J^|sgn[x 2 (t) - Xj(t)J + sgnjx 2 (t) - Xg(t) jdt d£ (116b) 


x 3 (*) = x 3 ( t= 0) + v 3 (t=0)t + 47rGm J^|sgn|x 3 (t) - XjftjJ + sgnjxg(t) - x 2 (t) jdt d? (116c) 


If a coordinate system is chosen such that 

Xj(t=0) + x 2 (t=0) + x 3 (t=0) = 0 

and 

v^(t=0) + v 2 (t=0) + Vg(t=0) = 0 

then 

Xj(t) + Xg(t) + Xg(t) = 0 

Thus, whenever the coordinate x 2 equals the coordinate x 3 , a straight line 

Xj = 2x 3 = 2x 2 


(117a) 

(117b) 

(118) 


(119) 


is obtained and a near-periodic behavior of the system is to be expected. Figure 45 
shows the results for a four-particle system. Again, the variation of the time-averaged 
velocity distribution and density indicates that there is near-oscillatory behavior of the 
system. This near-oscillatory behavior is also observed in plots of the kinetic-energy 
correlation function 


J 
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C(t) = (\T(t) - <T>][T(t+r) - <T>J> 

and also by plotting directly the kinetic energy of the system as a function of time. Fig- 
ure 46 shows the energy correlation function as a function of r for a six-particle sys- 
tem. The figure shows that the energy correlation decays only very slowly with r, an 
indication of near-oscillatory behavior. The fluctuations of the kinetic energy of the 
systems investigated are very violent and show no decrease with time. Even though 
these fluctuations show no decrease with time, the computer results show that a steady 
distribution can be reached while the fluctuations in the kinetic energy remain extremely 
violent. The time- averaged velocity distribution and density for the six-particle system 
are shown in figure 47. The experimental points are near the Maxwellian distribution 
for both curves. Figure 48 shows that for a 10-particle system the experimental velocity 
distribution and density points reproduce the theoretical Maxwellian curves. 

The results presented so far were obtained by time averaging the system over times 
of the order of N^t c , and they indicate that the time- averaged distribution is very nearly 
Maxwellian. However, to show relaxation to a Maxwellian distribution, a system should 
be prepared which is initially non- Maxwellian and then the relaxation to a Maxwellian 
distribution observed. This procedure was followed and the results for a 40-particle 
system are shown in figure 49. The solid line in the upper plot corresponds to a 
Maxwellian velocity distribution. The circles represent the initial velocity distribution 
of the system obtained by time averaging over 5Nr c . After the system was advanced in 
time for 2N^ characteristic periods, it was found that the experimental velocity distri- 
bution was still nearly the same as the initial distribution. Similar results are obtained 
for the density shown in the lower plot of figure 49. These results would indicate that no 
relaxation to a Maxwellian distribution occurs to order N^t c , but to be conclusive the 
calculations should be extended to 5N^ or 10N^ characteristic periods. The calcula- 
tions were terminated after 2N^ periods because of computer time limitations. Fig- 
ure 50 shows the time variation of the kinetic energy for the 40-particle system. The 
corresponding energy correlation function is shown in figure 51. The correlation of the 
kinetic energy appears to be very persistent and fluctuations seem to be superimposed on 
fluctuations. For large r the results resemble an amplitude modulated wave. The two 
cases shown in figure 51 were obtained by averaging over and 2N^ characteristic 

periods. It can be seen that if the average is taken over a smaller time interval, then 
modulation of the basic frequency for large t is much more pronounced whereas for 
small r the results do not change. The results presented in figure 51 indicate that for 
small r there exists Landau damping but some smaller correlation persists for a long 
time. 
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It is also of interest to determine whether the fluctuation of the kinetic energy of the 
one -dimensional system is due to thermal fluctuations. This has been done by plottings 
the relative fluctuation of the kinetic energy 


_ \lcj 0 ) 

<T > 


( 121 ) 


as a function of the total number of particles in a system N. The results are shown in 
figure 52. The solid line corresponds to i± = l/]lN and as can be seen, most of the 
experimental points are near the solid line. Thus the magnitude of the fluctuations of 
the kinetic energy is of the same order as would be expected from thermal fluctuations. 
The same result has been obtained by M. Henon with his model of concentric spherical 
mass shells. 


The rate at which the system will achieve equipartition of energy among stars of dif- 
ferent mass has also been investigated for systems containing sheets of mass m and of 
mass 4m. It was found to be convenient to present the results by means of the quantities 

= ( <T>r )’ 1 J q T t i dt 




( 122 ) 


: 2 .(<T>1)-‘J o T T 2 


dt 


where Tj is the total kinetic energy of the heavy stars and T 2 is the total kinetic 
energy of the light stars. The upper plot in figure 53 shows the results for a 10-particle 
system containing 5 heavy and 5 light stars. The time-averaged kinetic energies for the 
heavy and light stars approach each other. This result indicates that equipartition of 
energy occurs. Similar results are shown in the lower plot in figure 53 for a system 
containing 5 heavy and 15 light stars. The time -averaged kinetic energies of the heavy 
and light stars are now in the ratio 5/15, another indication that equipartition of energy 
occurs to order N^r c . 


CONCLUDING REMARKS 

The equilibrium properties of one-dimensional self-gravitating systems are inves- 
tigated analytically. One- dimensional computer models are used to perform experiments 
tracing the evolution of stellar systems. 

For the one-dimensional stellar system it was shown that for an initial distribution 
of the waterbag type with a given initial area in phase space the equilibrium state is a 
minimum-energy configuration. Numerical experiments with large numbers of particles 
illustrate these results. Thus, starting from any nonequilibrium state, the stationary 
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state described by the time -independent Vlasov equation can never be completely reached 
because the excess initial energy causes spiral arms to develop in phase space. These 
arms are required to accommodate the excess energy, and they wind around the main 
system and become longer and thinner as time progresses. After a sufficiently long time 
it becomes increasingly more difficult to detect changes in the distribution function since 
successive turns of the arms approach each other. It may then be possible to obtain an 
average distribution function by averaging over the gravitational potential. Such an anal- 
ysis has been performed by Lynden-Bell (Monthly Notices of the Royal Astronomical 
Society, vol. 136, no. 1, 1967). Lynden-Bell suggests that the encounterless relaxation 
of an unsteady system will lead to an equilibrium related to Fermi-Dirac statistics. The 
results of the present investigation of encounterless relaxation with a one-dimensional 
model appear to confirm the Lynden-Bell theory. Systems with an initial waterbag dis- 
tribution near equilibrium were found to approach a Fermi-Dirac distribution closely 
with the exception of a high-energy concentration of particles. For initial distributions 
far from equilibrium the system approaches a Maxwellian distribution as predicted by 
the, Lynden-Bell theory. This result should not be too surprising since relaxing condi- 
tions which lead to the distribution given by Lynden-Bell are expected to persist only 
over regions where the system is dynamically unstable. Computer calculations show that 
for initial distributions other than the waterbag type the results are nevertheless very 
similar to those for the waterbag distribution. An extension of the minimum- energy 
property showed that all one-dimensional stationary distributions that are always 
decreasing in going outward from the center of the system are stable. Numerical experi- 
ments were performed for two systems which did not satisfy this criterion and these two 
systems were found to be unstable. 

The waterbag model is found to have the interesting property that for initial energies 
near the equilibrium energy the system approaches its equilibrium configuration very 
closely. For initial energies far from the equilibrium energy the system is unstable and 
breaks up into smaller clusters. 

The one-dimensional model is also of interest as an approximation to the distribu- 
tion of velocity and mass normal to the galactic plane of a greatly flattened galactic sys- 
tem. For the two distributions investigated in this study the theoretical gravitational 
field E(x) normal to the galactic plane was compared with the results deduced by Oort 
from observations. The variation of E(x) for the Maxwellian distribution was found to 
be in good agreement with the experimentally obtained variation. The difference for 
small values of the position coordinate x is probably due to the central core of the 
galaxy. As expected, for the waterbag distribution the force E(x) differs markedly 
from the experimentally obtained variation for large x. 
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From the results obtained in the study of thermalization effects in a one- dimensional- 
model for a stellar system, relaxation to a Maxwellian distribution appears to take place 
to times of the order of N^r c where N is the total number of particles in the system 
and t c is the characteristic period. These results were obtained by taking a time 
average for systems containing small numbers of stars. It was also found that equiparti- 
tion of energy among sheets of different mass occurs to order N^t c . 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., March 25, 1968, 

129-02-01-01-23. 
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Figure 9.- Time development of kinetic energy for two equivalent systems with different graininess. 
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Figure 13.- Approach to equilibrium in phase space for system of 2000 mass sheets with ratio of initial to equilibrium energy equal to 1.33. 



Figure 13.- Concluded. 


Figure 14.- Variation of period of oscillation of mass sheets with energy 
for equilibrium waterbag model. Time t 0 has been arbitrarily chosen 
































Figure 30.- Time development of density 



X 


unstable system, 




r 



.30 . 
.25 . 
.20 . 


f(v) 


.101 


.05 L 


oL_ 

-5 


t * 142. 5r 




Figure 31.- Time 


t = 117.5r 

c ,o 




t = 197. 5 t 

C,0 






for unstable system, 




















3xl0 3 

2 - 



x X 

Figure 43.- illustration of time development of unstable two-contour system with v = 0.25. 
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Figure 4#.- Time-averaged velocity distribution and density for 10-particle system. 
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Figure 49.- Time-averaged density and velocity distribution for 40-particle system. 
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